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We investigate the dynamics of r modes at amplitudes in the nonlinear regime for rapidly but 
uniformly rotating neutron stars with a polytropic equation of state. For this, we perform three- 
dimensional relativistic hydrodynamical simulations, making the simplifying assumption of a fixed 
spacetime. We find that for initial dimensionless amplitudes around three, r modes decay on 
timescales around ten oscillation periods, while at amplitudes of order unity, they are stable over 
the evolution timescale. Together with the decay, a strong difi'erential rotation develops, conserving 
the total angular momentum, with angular velocities in the range 0.5. . . 1.2 of the initial one. By 
comparing two models, we found that increasing rotation slows down the r-mode decay. We present 
r-mode eigenfunctions and frequencies, and compare them to known analytic results for slowly ro- 
tating Newtonian stars. As a diagnostic tool, we discuss conserved energy and angular momentum 
for the case of a fixed axisymmetric background metric and introduce a measure for the energy of 
non-axisymmetric fluid oscillation modes. 
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I. INTRODUCTION 

The r modes of neutron stars are a subclass of inertial 
modes, i.e. modes where the Coriolis force is the main 
restoring force, which have purely axial parity in the slow 
rotation limit. In the limit of slowly and rigidly rotating 
Newtonian stars, the r-mode eigenfunction and frequency 
is analytically known. Its frequency, as for all inertial 
modes, is proportional to the angular velocity O of the 
star. The r modes are current dominated, that is the 
amplitude of the density perturbations is smaller than the 
velocity perturbations by a factor Q.. Because of this and 
also due to the lower frequencies, they are weak emitters 
of gravitational waves compared to modes of same energy 
for which pressure is the dominant restoring force. 

The extension of the r-mode solution to the case of 
rapid rotation and/or General Relativity (GR) is an on- 
going field of research. In the slow rotation approxi- 
mation, GR equations have first been derived by [J [5] . 
In [Hj it was shown that they are valid only for non- 
barotropic equations of state (EOS), and that the equa- 
tions for barotropic stars are qualitatively different. Fur- 
ther corrections to the equations were pointed out by [3] . 
The equations in [H H] for non-barotropic stars contain 
singular points, whose interpretation is still under de- 
bate [3JISHI2]. In particular, |5j provides arguments that 
the slow rotation approximation breaks down near the 
singular points, and including higher order terms would 
lead to valid physical solutions. However, there is still 
no mathematical proof that regular r modes of non- 
barotropic stars in GR exist under all circumstances. For 
the barotropic case, i.e. an EOS where the pressure is a 
function of density alone, GR counterparts of the New- 
tonian r mode have been found by [Sj. Those solutions 
are hybrid modes containing both polar and axial parity 
perturbations. 

The rapidly rigidly rotating case has been investigated 



in GR only with the simplifying assumption of a fixed 
spacetime (relativistic Cowling approximation). In (13j . 
the two-dimensional partial differential equations govern- 
ing the eigenfunctions for the case of a barotropic EOS 
have been solved numerically by using finite differences. 
It was found that the ratio of oscillation and rotation fre- 
quencies is decreasing with increasing ratio of rotational 
to binding energy. In [14] . the r modes and pressure 
modes of rapidly but rigidly rotating polytropic stars are 
investigated using time-evolution of the linearized equa- 
tions for a prescribed ^-dependency. The method is ex- 
tended to include the non-barotropic case in [15. , and 
differential rotation in [T5]. 

The r mode is of astrophysical interest since in GR it 
is generally subject to unstable growth, emitting gravita- 
tional radiation, as shown by [T71 [TB] . The reason is that 
the r-mode wave pattern is generally counter-rotating 
in the co-rotating frame, but co-rotating in the inertial 
frame, which is the condition for the CFS mechanism dis- 
covered by [ini [20] to be operational. Due to the CFS 
instability, r modes are a possible source for gravitational 
wave astronomy, but might also explain the limitation of 
observed neutron star rotation rates. It is thus important 
to know at which amplitudes the r-mode instability satu- 
rates. There are several effects which might suppress the 
instability or limit the amplitude to small values, includ- 
ing nonlinear couplings |21H23j . winding up of magnetic 
field lines by differential rotation associated with the r 
mode [24ll26j . and bulk viscosity [27] enhanced by the 
presence of hyperons. The studies mentioned so far ap- 
ply to newborn neutron stars, which are still hot enough 
to prevent superfluidity, which further complicates the 
picture, see [28H30] . 

Assuming the r mode has an instability window, i.e. 
a range of temperature and rotation rate where the CFS 
instability is not suppressed, it would be important to 
know the behavior at very high amplitudes, typically ex- 
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pressed as a dimensionless quantity proportional to the 
ratio of velocity perturbation to rotational velocity at 
the equator. Using numerical simulations of uniformly 
rotating stars in the Newtonian framework, it was found 
by that formation of shocks and breaking of surface 
waves set in at dimensionless amplitudes around three. 
However, the r modes in these simulations have been con- 
tinuously excited by a gravitational back-reaction force 
artificially increased to huge values. In a similar sim- 
ulation [231 E2], but switching off the artificial driving 
force after amplitudes in the range 1.6... 2. 2 had been 
reached, a catastrophic decay of the r modes was found. 
This decay was attributed to mode coupling effects, due 
to the observed growth of secondary modes. In full GR, 
a numerical simulation of an r mode with unit amplitude 
was performed by 331 , but no decay apart from numer- 
ical damping was found on the evolution timescale of 26 
oscillation periods. 

The r modes are also linked to differential rotation in 
many ways. For the case of rapidly and differentially ro- 
tating barotropic Newtonian stars, [M* found that the r- 
mode eigenvalue problem for a polytropic Newtonian star 
becomes difficult to solve numerically for degrees of dif- 
ferential rotation large enough to cause corotation points 
of the wave pattern. On the other hand, [35 found reg- 
ular solutions for a Newtonian, incompressible thin shell 
model with a similar amount of differential rotation. The 
existence of regular r-mode solutions for differentially ro- 
tating stars is thus likely, but remains unproven. For uni- 
formly rotating stars, the presence of r modes implies a 
certain amount of differential rotation, which is of sec- 
ond order in the mode amplitude. In [211 US], this effect 
was estimated using only the first order eigenfunctions. 
The importance for the evolution of a magnetic field via 
winding up of field lines was demonstrated in [24] |26] . 
Complementary to this, it was shown by |36| that the sec- 
ond order corrections to the eigenfunction of the velocity 
field necessarily contain a differentially rotating part, at 
least for slowly and rigidly rotating Newtonian stars with 
barotropic EOS. As argued in [371135]) the resulting sec- 
ond order contributions to the angular momentum have 
to be considered for the time evolution of the r-mode 
instability. A similar, but distinct effect concerns the in- 
teraction with gravitational radiation. For a toy model 
consisting of a spherical shell, it was shown by [33] that 
the gravitational back-reaction directly induces differen- 
tial rotation instead of causing only a uniform spindown. 
Further, the aforementioned studies [31] [32] of nonlinear 
r-mode decay also observed the formation of significant 
differential rotation during the decay. These studies sup- 
port the view that the CFS instability of the r mode will 
cause a certain amount of differential rotation, even if 
the r-mode amplitude saturates at small values. 

In the present work, we extent the previous studies in 
mainly two directions. First, we extract r-mode eigen- 
functions and frequencies for rapidly (but rigidly) rotat- 
ing relativistic stars under the simplifying Cowling ap- 
proximation. We are using a fundamentally different ap- 



proach than [T31 [U] , and also discuss the properties of 
the eigenfunctions in detail, in particular their energy 
and estimated gravitational radiation. Second, we inves- 
tigate the decay of high amplitude r modes and the for- 
mation of differential rotation already found in [3T] [32] . 
We are however using the relativistic Cowling approx- 
imation instead of Newtonian gravity for the evolution 
and the exact (linear) r-mode eigenfunctions as initial 
perturbation, instead of breeding r modes by using ar- 
tificial back-reaction forces. This will shed light on the 
robustness of those effects. 

II. ANALYTIC TOOLS 

In this section, we review the analytic tools used for 
this work. Throughout the article, we use the following 
notation: 



p 


= Rest frame rest mass density 


p 


= Pressure 


e 


= Specific internal energy 




P 


h 


= ! + €+- 




P 




= Fluid 4-velocity 


v' 


= Fluid 3-velocity 


W 


= Fluid Lorentz factor 




= 4- metric of signature (—, + ,+,-!-) 


9ij 


= 3-metric 


a 


= Lapse function 




= Shift vector 




= Angular velocity of the star 



All equations assume geometric units G — c = 1. Greek 
indices run from ... 3, indices i,j,k,l from 1 ... 3. When 
working in cylindrical coordinates, generally denoted by 
("i?, z, 0), indices a, 6, c run from 1 ... 2, excluding 0. 

A. General properties of eigenfunctions 

Although there is no analytic solution for the oscil- 
lation eigenfunctions of rapidly rotating stars, one can 
derive some general properties. The following is valid for 
a fixed spacetime, rigid rotation, and a barotropic EOS. 

For any global oscillation mode, the perturbation of a 
fiuid quantity X can be written as 

SX {x, t) = [x z) e^C^t+^-^+Sx)^ ^ 

where A is a dimensionless amplitude, w the real-valued 
oscillation frequency, and Ox is a constant phase shift. 
X is the (suitably normalized) real-valued eigenfunction. 

The set of variables we use to completely specify the 
oscillation is {e, w*}. The relative phase shifts are given 
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by 9^ = 9y4, = 0, 0i,t> = Oyz = |. For a derivation, 
see [to]. 

Due to the equatorial symmetry of the unperturbed 
models, the eigenfunctions also have a well defined z- 
parity, with the relations 



P,[Se] = P^v"] = ~P^[Sv^] - P.M. 



B. Newtonian r mode 



(2) 



For slowly rotating stars in Newtonian theory, there 
exist analytic solutions for the r modes. The velocity 
eigenfunctions are given by 



where 



Yil = {l{l + l))--rxVYir. 



(3) 



(4) 



are the pure-spin vector harmonics of magnetic type (see 
e.g. III]). The density perturbation is proportional to 
Q^Yz+i^i, see [42], and the frequencies oji in the inertial 
and Wc in the corotating frame are 



i + i 



LOc 



2n 

ITT' 



(5) 



see [43] . The negative sign of uji means prograde motion 
of the wave patterns. 

To measure the r-mode amplitude during a simula- 
tion, it is common to use a scalar product with the 
magnetic vector harmonic times some radial weighting 
function. This makes sense because the scalar product 
with other oscillation modes vanishes in the slow rota- 
tion limit. A convenient choice to measure the r-mode 
amplitude is given by the magnetic current multipole mo- 
ment Jii , which is also used to estimate the gravitational 
wave strain. For the Newtonian r mode, we find 



5Ju - j pr'5v-Yif*d^x, 



yB 



Hence 



A = 



\SJu\ 



yB 



d3a 



(6) 
(7) 

(8) 



In this work, we use the above formula to define the di- 
mensionless r-mode amplitude also for the rotating rel- 
ativistic case, setting R to the circumferential equato- 
rial radius, and evaluating the denominator for the back- 
ground model. Note this is not a covariant measure, for 
exact comparison to other works one should use invari- 
ants like total energy or maximum velocity. With in- 
creasing rotation rate, one can expect that the presence 



of other modes (with the same m) starts contributing to 
the above measure as well. This should not be a problem 
as long as the r mode is the dominant one. Differential 
rotation and other axisymmetric perturbations do not 
contribute to Ju. 



C. Evolution equations 

We evolve the general relativistic hydrodynamic equa- 
tions for an ideal fiuid, which in covariant form read 



(pu^) = 0. 



(9) 
(10) 



The stress energy tensor T'"' of an ideal fluid is given by 

T^""" ^ phu^'u" + Pg^"" . (11) 

For numerical evolution, a 3-f 1-split is applied to obtain 
a first order system of hyperbolic evolution equations in 
conservation form with source terms 

d^q^~d,nq,x') + s{q,x'), (12) 
q={D,T,S,), (13) 

with the evolved hydrodynamic variables given by 

D = ^Wp, (14) 
T = ^{W^ph-P-Wp), (15) 
= ^W^phv,. (16) 

In flat Minkowski space (and using standard coordi- 
nates), _D, r, and S** reduce to mass density, energy 
density not including rest mass, and linear momentum 
density. The flux terms p = {f}^, f^, fg.) are given by 



Pd = w'D, 
n = + a^v'P, 
Ps^=w^S,+a^P5], 



(17) 
(18) 
(19) 



where = aw' — /3* is the advection speed relative to 
the coordinate system. The source terms can be writ- 
ten in many ways, the formulation we are using is dis- 
cussed in [33]. Finally, the evolution equations need to 
be completed by an equation of state (EOS) of the form 
P — P{p, e) to compute the pressure. 



D. Conserved quantities 

Making the assumption of a fixed axisymmetric space- 
time not only simplifies the numerical evolution, it also 
implies the existence of conserved fluxes. The station- 
arity of the metric leads to a conserved energy density. 
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while the the axisymmetry of the metric leads to a con- 
served angular momentum. Conserved mass, energy, and 
angular momentum are given by 



M 



Dd^a 



E = / Ud?x, J = / Ld^x, (20) 



where D is defined by Eq. ( 14 ) and, using coordinates for 



which 9t,90 are Killing vectors. 



U = ^T^„n'' -D 

aP' 
Wp, 

^T^un" = DWhv^. 



= D I Wha - 1 

L — Stj, 



(21) 
(22) 
(23) 



Above, denotes the unit normal to the surfaces of 
constant coordinate time. The conservation of M, i5, J 
can then be easily derived from the evolution equations 
(12 1, without using the Einstein equations. For an arti- 



ficially fixed spacetime, conservation thus holds even for 
non-axisymmetric, non-stationary fluid flows. The only 
requirement is that the fixed background metric is sta- 
tionary and axisymmetric. 

Note that E depends on the gauge quantities a and 
p. Since we require that dt is a Killing vector, the only 
freedom of choice for j3 is given by translations along 
space-like Killing vectors. The resulting change amounts 
to adding multiples of the conserved quantities associ- 
ated with those, e.g. angular momentum for a rotation. 
The only freedom of choice for a is multiplication by a 
constant, which results in replacing i? by a linear com- 
bination with M. For the rest of this article, we assume 
asymptotic fiatness and define E as the energy computed 
for a coordinate frame where and a — )• 1 at spa- 

tial infinity. When making this choice, C/ — for an in- 
finitely diluted fluid element at rest at infinity, and hence 
a system becomes unbound for i? > 0. 

Also note that the numerically evolved energy density 
r is not the same as the conserved energy density U . In 
the Newtonian limit for example, for a fiuid moving in 
a constant external gravitational field, t does not con- 
tain potential energy, in contrast to U . The reason for 
not using U in numerical schemes is that the correspond- 
ing evolution equation involves the time derivative of the 
gauge quantity a in the fully relativistic case. 

We stress that only M is still conserved in full GR, and 
E is not even conserved in Newtonian physics when al- 
lowing the gravitational field to change. Further, E + M 
is different from the ADM and Komar mass. Interest- 
ingly, J equals the Komar angular momentum, as long 
as there are no singularities to consider. 



E. Oscillation energy 

In the following, we introduce a notation for the en- 
ergy of oscillation modes in the Cowling approximation, 
which is useful as a diagnostic tool in numerical studies 



performed in a fixed spacetime, but might also serve as 
an order of magnitude estimate in the general case. 

For this we expand the conserved energy E defined 
in the previous section around a stationary background 
model. For non-axisymmetric modes, it is easy to show 
that the linear terms cancel when integrating along the 
(/i-direction, so we need to expand the energy to second 
order in the amplitude. Since we want an expression that 
can be computed without knowledge of the second order 
corrections to the eigenfunction, we still assume that the 
perturbation itself scales linearly. However, this leads 
to conceptional problems. First, there is a freedom of 
choice concerning the set of variables used to completely 
specify the system. If those variables are perturbed lin- 
early, the resulting energy perturbation at second order 
depends on this choice. Second, we found that in general 
the conserved mass changes as well when evaluated to 
second order. This implies that the state around which 
the model oscillates is not exactly the same as the un- 
perturbed state. Unfortunately, we are interested in the 
energy difference from the stationary state reached after 
all the oscillation energy is dissipated, e.g. due to numer- 
ical damping, while conserving total mass and angular 
momentum profile. For one of our models, we computed 
mass and energy change when perturbing p, with the 
r-mode eigenfunction. By assuming that the mass is cre- 
ated with the average binding energy E/M oi the back- 
ground model, we estimated the ambiguity in defining E 
to be around 50 %. This effect could be tracked down 
to the fact that not only the kinetic energy density, but 
also the conserved mass density depend on the velocity 
via the Lorentz factor, in contrast to the Newtonian case. 

To cure these problems, we define the oscillation en- 
ergy as the perturbation of the conserved energy when 
perturbing the variables £',L,ti" linearly. This way, to- 
tal mass and angular momentum are exactly conserved 
for non-axisymmetric perturbations. To our knowledge, 
the following has not been discussed elsewhere. We will 
assume a cylindrical coordinate system adapted to the 
symmetries and with /S" — ga,f, = 0, and only consider 
background models with v"^ — Q. For a non-axisymmetric 
mode described by Eq. ([T]), we define the mode energy 
as 



E 



5D^ 



+ 



SDSL 



2dL^ 



6L' 



dDdL 

1 d^U 



2 91)2 



dDdL 



DL 



1 



2 dv'^dv'' 
^ 2 



Sv'^Sv'' d^x 



1 



2 dv'^dv'' 



(24) 



(25) 



v^'v'' I dddz. 



The quantities in Eq. (24) are defined with respect to 



Cartesian coordinates, while Eq. (25) is valid in cylin- 
drical coordinates. Terms with SDSv"' or SLSv"^ do 
not contribute since they have the angular dependency 
sm{m(j)+ijt) cos{m(f)+ujt). Also the corresponding second 



derivatives are zero. This energy depends on the normal- 
ization of the eigenfunctions, which has to be specified. 
The eigenfunctions of the conserved variables are given 
by 

D = {p + pW'^v^v^) , (26) 

L = ^W^h [(1 + cl) v^p + p {2W^ - 1) g^^v'f'] . (27) 

To compute the second derivatives, we first define 

U{W, p, L) = -P'I'L + ^ {Wp (Wha - 1) - aP) (28) 
= U{D,L,v''). (29) 



rive at 



We then write 



5L2 



dDdL 



(PlJ_ /dW 
dm \~dD 



&HJ_ 
dp^ 



dD 



d^U fdwy 

:)2fr / o N 2 



9p2 



dp 
dL 



d^U dW dW 



+ 



dU d^W dU d^p 
dW dD^ ^ ~d^~dlfl 

d^U dp dW 
dWdpdoW)' 

d^tJ dpdW 
dWdpdLl)L 

dU d^W dU d^p 
' dW~dlJ ^ 'd^dL^' 
d'^U dp dp 



dW^ dD dL dp^ dL dD 
d^tJ f dp_dW_ dW dp 
dL~dD ~^ 'dLdD 
dU d'^p 



dWdp 
dU d^W 



dp dDdL' 
dtJ d^p 

= 1 — , 

QyaQyb Q^^ QyaQyb Qp QyaQyb ' 



d^u 



dW dDdL 
dU d^W 



(30) 



(31) 



(32) 



(33) 



where we included only nonzero terms. From Eq. (28 1, 
we compute 



dn'l' 




yj J-J 


y/^VV pj 




1 

i 


~dL " 


^W^phf 


dp 


2W^ - 1 


dD 




dp 




dL 





where 



/ = W'vl (1 - cl) + 



It follows that 



dW W^v'^ 



dD ^jpf 
dW _ W 
dL yffphf'^ 



Using Eq. (39)-(45l, we obtain 



dl 
dD 



dL-'^'^-^hf 



V7/ 



m 

I 

-(1- 

p 



2d 2 



Finally, we find 



dU 
dW 
dU 

'dp 
(PU_ 
dW^ 
d^U 



dp"^ 

dm 

dWdp 



^p {2ahW ~ 1) , 

{ahW (1 + v^cl) - 1 
2^aph, 

^ {2ahW {I + cl) ^ I) . 



2 d 2 



(34) 
(35) 
(36) 
(37) 

(38) 



The functions W{D, L,v°-), p{D, L,v'^) cannot be ex- 
pressed in closed analytic form. However, we only need 
the derivatives. By computing the derivatives of the con- 
served variables with respect to the primitives and then 
inverting the resulting linear system of equations, we ar- 



dD^ 



d^W 



dLdD 



dv°'dv^ 



YdD-'"^ ""^dD 



1 + ci 



_ d_ . 
" dp"'' 



W 



dL^ ^phf 
V4, df 



dv'l' 
~dL 



9p 
dD 



dp 



m 

Vipf 



f dL p dL 



W 



/iphf 



dv* 



v^df_ 2\ dp 

' f dD p ^ "^''^^ dD 



T 



-9<p<pgab, 
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9V 



dLOD 
JPp_ 

dyaQyb 



9<t><t> 

- {2W^ - 1 
1 



W {2W^ + 1) 



dv't' 



1 

7w 



_j_di 

Wf dD 

V4> df 
f dL 

f dD 



- 94>4> 



dv't' 



p dL 

4 dp 
p dD 



f 



-pW^gab- 



(52) 

(53) 
(54) 
(55) 



First and second order derivatives have of course been 
computed without assuming v°' — 0, but the results given 
here have been evaluated at v°' — for simplicity. 



We stress that the energy defined by Eq. (24 1 is only 



an estimate for the energy of a finite amplitude oscilla- 
tion. It was shown in 05] that for r modes of Newto- 
nian stars, differential rotation is an unavoidable feature. 
Most likely, such axisymmetric terms would contribute 
to the first order expansion of the energy E and hence 
constitute a second-order contribution in total, like the 
terms considered in our definition. Without knowledge of 
the second-order perturbation, Eq. (24 1 is probably the 
best one can do. 



To evolve the above system numerically, we use the 
PIZZA code first described in [44j . It is based on an HRSC 
(high resolution shock capturing) scheme, which was 
optimized for quasi-stationary simulations. As shown 
in [33], the code is able to evolve a stationary star with 
high accuracy, in particular the rotation profile. As for 
all such codes, special care has to be taken to treat the 
stellar surface. Instead of using an artificial atmosphere, 
we apply the scheme used in [31]. The advantages are 
that the mass is conserved to machine precision and that 
the amplitude of oscillations excited by numerical noise 
is negligible for the applications presented here. Still, the 
surface treatment is the main source of numerical damp- 
ing. 

We use uniform three-dimensional Cartesian grids in 
co-rotating coordinates. The outer boundaries are placed 
far enough from the surface to prevent matter from leav- 
ing the computational domain. Thus, the total mass M 
is conserved to machine precision. The total angular mo- 
mentum on the other hand is subject to discretization 
errors, since L is not an evolved variable in Cartesian 
coordinates. 



III. NUMERICAL METHOD 



B. Eigenfunction extraction 



In the following we briefly describe our numerical 
methods. For readers not familiar with general relativis- 
tic hydrodynamics, we recommend the review [3S]. 



A. Time evolution 

We evolve the 3-1-1 split hydrodynamic evolution equa- 
tions in flux-conservative form (12). However, we use a 
zero-temperature EOS of the form P — P{p). As a conse- 
quence, the system becomes overdetermined. Therefore, 



we do not evolve the energy density r, which becomes 
redundant, but compute it from mass and momentum 
densities. For details, see |33j- We stress that this ap- 
proach is only self-consistent for adiabatic evolution, and 
hence our simulations are only correct in the absence of 
shock formation. Discontinuities cannot produce shock 
heating, instead they lead to a violation of energy con- 
servation. A sharp decrease of E is therefore an indicator 
for shock formation, see for examples. 

In the absence of shocks however, the evolution be- 
comes more accurate, since there is no error in the evolu- 
tion of the specific entropy. One particular error that is 
avoided this way is the formation of large scale, low ve- 
locity convective movements driven by entropy gradients 
caused by numerical errors. Such vortices could easily be 
confused with genuine nonlinear effects in simulations of 
high amplitude r mode oscillations. 



In order to extract eigenf unctions numerically, we use 
the mode recycling method in the form described in j46| . 
In short, the star is perturbed using some trial perturba- 
tion, and then evolved in time numerically. The frequen- 
cies of the oscillations excited by the perturbation are 
determined using Fourier analysis in time at some sam- 
ple point. By selecting one frequency and computing 
the Fourier integrals at this frequency for any point in the 
star, one obtains a first estimate of the complex eigen- 
function, which is usually still contaminated by other os- 
cillations due to the finite evolution time. To obtain the 
two-dimensional eigenfunctions, we divide the numeri- 
cal complex eigenfunction by e"""^ and average over (p. 
This removes contributions with different 0-dependency. 
Next, we remove contributions with the wrong z-parity 
by (anti-)symmetrizing. Finally, we compute and factor 
out the average complex phase (see 0^ for details), while 
using the phase variance as an error measure. The result 
is used as the initial perturbation in a new simulation, 
and the whole process is repeated until oscillation modes 
other than the desired one are reasonably suppressed. 

Obviously, the method is only effective if the initial 
trial perturbation significantly excites the desired mode. 
To extract the r-mode eigenfunction, we choose the New- 
tonian eigenfunction for a slowly rotating star, given by 
Eq. ([3|, which is close enough to the actual eigenfunc- 
tions of the models. 
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Name 


Mb/Mq 




_Rc/km 


ar 


fie 


/3c 


Model 


./c/Hz 


/./Hz 




E/E 


MB85 


1.6194 


590.90 


15.384 


0.85 


-0.2094 


0.02406 


MB85 


361.4 


820.4 


0.917 


3.033-10-'' 


MB70 


1.7555 


792.10 


17.268 


0.7 


-0.2159 


0.04864 


MB70 


544.8 


1039.4 


1.032 


5.093-10-'' 



TABLE I: Details of the stellar models. Mb is the total 
rest mass, Fb the rotation rate as observed from infinity, Rc 
the equatorial circumferential radius, the ratio of polar to 
equatorial coordinate radius, = E/Mb, Pc = \^J/{2E)\, 
where E and J are energy and angular momentum defined by 
Eq. pOl. 



IV. STELLAR MODELS 



TABLE IL Properties of the I = m = 2 r-mode. fc is the fre- 
quency with respect to the co-rotating frame, fi the frequency 
observed from infinity in the inertial frame. The estimated 
numerical accuracy of fc is 1%, not including the unknown 
error due to the Cowling approximation. For comparison we 
compute the r-mode frequency for slowly rotating Newto- 
nian stars given by Eq. ([5|. -E is the r-mode energy, defined 
by Eq. ( |24[ ), normalized to an amplitude A = 1 as defined by 
Q. -E is the models' binding energy defined by Eq. (20 1. 



Eq. 



We investigate two different uniformly rotating 
stellar models with fixed central density pc — 
7.9053-10" kgm-3 and EOS, but different rotation 
rates. Their properties are summarized in Tab. |T] The 
EOS, which is also used during the evolution, is a poly- 
tropic EOS defined by 

m = P.{£)\ (56) 

with polytropic exponent T = 2 and the constant density 
scale pp = 6.1760 - 10^^ kg m-^. We note that polytropic 
stars are stable against convection and do not possess g 
modes, i.e. modes for which buoyancy is the restoring 
force. Model MB85 was computed using the rns code 
described in [17] , while model MB70 was computed using 
the code described in (351 133] - Both codes are accurate 
enough for our purposes, the only reason to use different 
codes is that the latter became our standard choice once 
it was available. 

These models are a crude approximation to real neu- 
tron stars. Their purpose is to get a basic qualitative 
understanding of the nonlinear r mode dynamics in the 
most simple case. It is likely that the inclusion of compo- 
sition gradients or differential rotation will lead to new 
effects. 

To excite nonlinear oscillations, we always use the ex- 
act eigenfunction obtained by mode recycling. However, 
in contrast to the linear regime, there is some arbitrari- 
ness involved regarding how to scale the eigenfunctions 
to large amplitudes. Ideally, we would like initial data re- 
sembling a mode naturally grown to high amplitudes, e.g. 
due to the CFS mechanism. Since this is not feasible, we 
simply scale the perturbation of the velocity and of the 
specific energy and recompute the other fluid quantities 
consistently after applying the perturbation. 

Note also that we use Eulerian perturbations. This 
has the drawback that the star is not deformed correctly 
close to the surface, in particular there is no perturbation 
outside the surface of the unperturbed star. For high am- 
plitudes, the initial data inevitably contains small shocks 
at the surface. In practice however, even high amplitude 
r modes do not induce large deformations and we did not 
notice the corresponding numerical artifacts. 



V. NUMERICAL RESULTS 

In the following, we present our results for the r mode 
with I = m = 2. Unless noted otherwise, our simula- 
tions use a uniform resolution of 50 points per equatorial 
stellar radius, which is a reasonable compromise between 
accuracy and computational cost. 

A. R-mode properties 

Using the methods in Sec. |IIIB[ we extracted eigen- 
functions and frequencies of the r mode for the models 
in Tab. |l] We also computed the energy of the modes 
defined in Sec. |IIE| The results are given in Tab. |Tlj For 
our models, the frequencies agree with the ones in the 
Newtonian slowly rotating case better than 10 %. The 
r-mode frequencies in the inertial frame found by |16j 
agree with our results better than 0.1%. Given that we 
use a completely independent code based on a different 
method, the good agreement validates the results. 

The mode energy is useful to quantify what amplitudes 
are large in the sense that strong deformations of the 
star occur. Naively, one should think that e.g. A = 3 is 
a huge amplitude because the velocity perturbations be- 
come comparable to the rotational velocity at the equa- 
tor. Looking at the energy however, we find that the 
mode energy at ^ = 3 for model MB85 is only a fraction 
3 - 10"^ of the binding energy E. For model MB85, the 
energy of the perturbation equals the stars' binding en- 
ergy aX A ^ 57. For the mode recycling process, we used 
amplitudes around A = 0.3, well inside the linear regime. 

The eigenfunctions are visualized in Figs. [l}l2j As one 
can see, the velocity perturbations are similar to those 
in the Newtonian slowly rotating case, in particular the 
radial component is small. Also the density perturba- 
tion is qualitatively the same as in the Newtonian case. 
To quantify the differences, we decompose the velocity 
perturbation into vector spherical harmonic functions 

oo / 

1=0 m=-l (57) 
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FIG. 1: Two-dimensional eigenfunctions e{'&,z) (left half) 
and t)"(i9, 2) (right half), belonging to the r mode of model 
MB85. 
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FIG. 2: Perturbations corresponding to the r mode of model 
MB85, as longitude/latitude plot at fixed coordinate radius 
r = Re/2, with being the equatorial coordinate radius. 
The velocity perturbations 5v'^,Sv^ are plotted on the left 
half, and the specific energy perturbation 5e on the right half. 
Note the patterns have a 180°-periodicity in ip. 

The results are plotted in Fig. [3j The dominant contribu- 
tion to the velocity is the magnetic type vector harmonic. 
However, there are also significant radial and polar con- 
tributions. The dominant component of the specific en- 
ergy (not shown in the plot) is the Z = 3, m = 2 spherical 
harmonic. 

Although higher order multipole moments are quite 
small in the inner regions of the star, they have signifi- 
cant amplitudes between the polar and equatorial radius. 
This does not imply that small scale structures appear 
close to the surface. The reason is just that the star is 
ellipsoidal and the spheres of constant coordinate radius 
start intersecting with the stellar surface. 

From the multipole decomposition of the numerically 
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FIG. 3: Decomposition of the r-mode velocity perturbation 
of model MB85 into vector spherical harmonics. The vertical 
lines mark the polar and equatorial radius. 



Model 


Aii/{W Mpc) 


te/s 


rj/s 


MB85 


1.42-10"^'* 


1.47- 


10" 


4.93-10* 


MB70 


3.56- 10-^"* 


1.64- 


10^ 


1.04-10* 



TABLE III: Gravitational radiation caused by r modes in the 
linear regime. The gravitational luminosity Wgm and angular 
momentum loss Jg^ axe given in terms of the timescales rj = 
J/Jgw and te ~ \E\/Wgjn, where E and J are the binding 
energy and angular momentum given in Tab.|l] The values are 
normalized to an amplitude A — I, with A22 ~ A, Wgw ~ A^, 

<J gw ^ . 

extracted eigenfunctions, we estimated the gravitational 
radiation caused by the r modes, using the multipole 
formulas for Newtonian sources from [41] . Luminosity, 
strain, and angular momentum loss are given in Tab. 
Unsurprisingly, we find that for the r modes only the I = 
m = 2 current multipole contributes significantly. The 
second largest contribution comes from the Z = 3, m = 2 
mass multipole, which for model MB85 is smaller by a 
factor \ A^2 /■^22 \ ~ 0.02. The higher multipole moments 
are completely negligible. 

B. Rotation profile 

One of the most noticeable features in our simulations 
is the development of strong differential rotation during 
the evolution of r modes with high initial amplitudes. 

To visualize the rotation profile alone without the con- 
tribution from the r-mode oscillation, we compute the 
(/)-averaged angular velocity. Of course, any axisymmet- 
ric oscillation would contribute to this measure as well. 
However, since the Fourier spectrum of the (/>-velocity at 
some sample points in the star did only show significant 
peaks at the r-mode frequency and at zero frequency, we 
can assume that any snapshot of the 0-averaged angular 
velocity is a good measure of the differential rotation. 
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FIG. 4: Development of difFerential rotation, for model MB85 
perturbed with an r mode of initial amplitude A = 3. Color 
coded is the (^-averaged angular velocity in units of the an- 
gular velocity of the unperturbed model. The surface of the 
unperturbed star is marked by a thin dotted line, while the 
surface at which the (^-averaged density falls below 0.005 times 
the central density is marked by a thin solid line. A small re- 
gion near the axis has been excluded to avoid problems when 
transforming from the Cartesian grid. 



Fig. [4] shows snapshots at two different times during 
and after the decay. During the decay, the angular ve- 
locity shows a two-dimensional structure. Strangely, the 
angular velocity near the poles temporarily increases to 
more than twice the initial value and then decreases 
again. At a later stage, the angular velocity converges 
to a simple profile depending roughly on the distance to 
the axis. In order to quantify the amount of angular mo- 
mentum redistribution, we define an average change of 
angular momentum by 



AJ 



L = 




(58) 
(59) 



where Lq is the value for the unperturbed model. Fig- 
ure [5] shows the time evolution of this measure as well as 
the local differential rotation at chosen positions. 

The profile shown in Fig. [4] is similar both in shape and 
magnitude to the differential rotation found by [31] for a 
Newtonian star. A cut in the equatorial plane is shown in 
Fig. [6) Near the axis, the rotation rate is slowed down by 
a factor of two, while the equatorial rate is increased by 
a factor 1.2. This also agrees well with the profile shown 
in [32] for the Newtonian case. It differs however strongly 
from the profile shown in [23]. On the other hand, the 
latter was not a (/)-average but a cut along the a;-axis. 
The magnitude of differential rotation is large enough to 
cause to a visible deformation of the stellar surface, as 
shown in Fig. |4| 

For the unperturbed models, our code is able to con- 
serve the rotation profile with errors several orders of 
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FIG. 5: Time evolution for the same setup as Fig. |4|of the 
0-averaged angular velocity at different sample points in the 
star. Also shown is the global average angular momentum 
change defined by Eq. (58 1. 
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FIG. 6; Snapshots at different times of the (/)-averaged angu- 
lar velocity in the equatorial plane versus coordinate radius, 
for the same setup as Fig. |4| The lines end where the (j)- 
averaged density falls below 0.005 times the central one. 



magnitude smaller than the observed differential rota- 
tion. The errors in the presence of differential rotation 
might be larger, because there is shear motion and be- 
cause the fluid is not co-rotating with the coordinates 
anymore, particularly near the stellar surface. In general, 
conserving the angular momentum is problematic with 
codes based on Cartesian grids. Therefore we monitor 
the conserved angular momentum J defined by Eq. ( 20 1 . 



For model MB85 with initial amplitude A — 3, the angu- 
lar momentum is decreasing more or less linearly, and the 
total loss at t = 20 ms is AJ/J = 0.0032. At the same 
time, we observe an average angular momentum change 
AJ/J = 0.26, much more than the total violation of an- 
gular momentum. It is also greater than the total and 
average change of angular momentum introduced by the 
initial perturbation (caused by second order terms and 
surface effects), which are of the order AJ/J — 0.018 
and AJ/J = 0.025. 

Although our computational resources did not permit 



10 



a 
<l 




3.5 



0.0 


























































































































' ' — 1 





10 15 20 25 30 35 40 45 

T/ms 



FIG. 7: Convergence of differential rotation. Shown is the (p- 
averaged change of angular velocity after 5 ms for an r mode 
of model MB85 with initial amplitude A = 3, at different grid 
resolutions. 



FIG. 8: Decay of the r-mode amplitude A defined by Eq. Q 
for model MB85 perturbed with different initial amplitudes. 
The two curves with initial amplitudes of ^4 = 2.48 and A = 
2.51 highlight the steep dependence of the decay time on the 
initial amplitude. 



a full convergence test, we evolved the first 5 ms with 
resolutions N = 37, 50, 75 points per stellar radius. For 
each we sampled the perturbation of 0-averaged angular 
velocity at the end of the simulation along the equatorial 
plane, as shown in Fig. [7] On average, low resolution 
seems to damp the differential rotation and not to cause 
it. To quantify the errors, we computed the i^-norms of 
the residuals and estimated the convergence order p ^ 
1.3, which at resolution A'^ = 50 implies an average error 
of AS7 around 10% of the maximum value. We note that 
the numerical evolution scheme is second order accurate, 
but the treatment of the surface is not. From Fig. [Tj one 
can see that the largest error of O is indeed found near 
the surface. 

We conclude that the differential rotation is caused 
by a redistribution of angular momentum, and not by 
numerical errors or contributions already present in the 
perturbed initial data. 



C. R-mode decay 

Simultaneously with the development of differential ro- 
tation, high amplitude r modes exhibit a rapid decay in 
amplitude, as will be shown in the following. 

As a measure for the decay, we use the dimensionless 
amplitude A defined by Eq. ( [8|. The results for all our 
simulations are shown in Figs. |8J and [9j All simulations 
were numerically stable (the shorter ones are exploratory 
simulations). 

As one can see, the decay depends crucially on the ini- 
tial amplitude. For initial amplitudes of order unity, the 
decay is comparable to the numerical decay, i.e. com- 
patible to no physical decay at all, while for amplitudes 
as high as 3, the amplitude decreases rapidly. Interest- 
ingly, for some initial amplitudes there is a plateau phase 
before the catastrophic decay. The length of this phase 
increases rapidly with decreasing initial amplitude. From 
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FIG. 9: Decay of the r-mode amplitude at different initial 
amplitudes, for model MB70. 



our simulations, it is unclear whether there is a critical 
amplitude for the onset of the decay, or if all r modes de- 
cay eventually. The decay we observe is qualitatively the 
same as the one reported in [32j , where a Newtonian star 
is evolved, exciting an r-mode by applying an artificially 
increased gravitational back-reaction that is switched off 
at given r-mode amplitude. 

The fact that the decay rate depends on the initial 
amplitude is a strong hint that the main cause is nei- 
ther wave breaking nor shock formation, since for those 
the damping strength typically depends on the instanta- 
neous amplitude, not on the initial one. By comparing 
Fig. [8] and Fig. [5] we can see that the increase of dif- 
ferential rotation is related to the r-mode amplitude. A 
hypothetical model which would explain the accelerating 
decay and saturating differential rotation is that differ- 
ential rotation causes r-mode decay and the presence of 
the r mode causes an increase of differential rotation. 

A comparison between models MB85 and MB70 shows 
that faster rotation stabilizes high amplitude modes. For 
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FIG. 10: Energy budget for a decaying r mode of model 
MB85 with initial amplitude ^ = 3. Shown is the loss of 
total conserved energy E defined by Eq. ( 20 1 , as well as the 
loss of energy in the r mode, computed by estimating the 
amplitude from the current multipole moment J22 and the 
mode energy from Eq. (24 1. For the points labeled "Phase 



ignored" , we took the absolute value of the complex integrand 
when computing J22. 



initial amplitudes of order unity, doing a comparison is 
difficult since the damping timescales becomes compara- 
ble to the numerical damping timescales for both models. 

To answer the question whether the energy of the r 
mode could be dissipated directly, e.g. via shock forma- 
tion or numerical dissipation, we now study the energy 
budget of the decay, using the tools from Sec. |IID| and 
Sec. |IIE| Fig. [To] shows the total conserved energy as 
well as an estimate for the energy in the r mode. The 
violation of energy conservation is obviously significant, 
but it is not sufficient to explain the amplitude decay, 
even if all the lost energy is taken from the r mode. We 
assume that the mode energy is at least partially trans- 
ferred to the differential rotation and deformation of the 
star described in Sec. IV Bl 

The only effects causing violation of energy conserva- 
tion are formation of shocks in conjunction with the cold 
EOS, surface effects like wave breaking, and numerical 
dissipation. To estimate the latter, we compare simula- 
tions of the first 5 ms at resolutions 37, 50, and 75 points 
per stellar radius. As can be seen in Fig.[TT] the violation 
of energy conservation is obviously caused for the most 
part by numerical errors, although some residual physical 
effect cannot be ruled out. However, the comparison also 
shows that the decay of the r mode is computed quite 
accurately. This implies that the loss of total energy, 
which does depend on resolution, is unrelated to the loss 
of r-mode energy. 

Still looking for the cause of the energy loss, we pro- 
duced a movie of the first 5 ms showing density and ve- 
locity perturbations in the coordinate planes. We did not 
notice any shock formation. The velocity field is domi- 
nated by the r mode. Right from the start, the density 
field is an overlay of many oscillation modes including 
the r mode and the quasi-radial mode. The presence of 
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FIG. 11: Convergence of r-mode amplitude (a) and conserved 
energy (b), for model MB85 perturbed with an r mode of 
initial amplitude A — 3. 



other modes is not surprising, on the contrary it would 
be strange if the simple linear scaling of the eigenfunction 
used to excite high amplitude r modes would not excite 
other modes as well. Further, the density perturbation 
of the r mode itself is weaker than for pressure modes 
of the same kinetic energy. We noticed an m = 4 defor- 
mation which near the surface becomes quite nonlinear, 
but not enough to cause wave breaking. However, the 
resulting numerical dissipation might explain part of the 
energy loss. In [31] . wave-breaking of the r mode itself 
was observed for a simulation of a Newtonian star after 
the r-mode amplitude peaked at A — 3.35. This does 
not contradict our results. Our maximum amplitude is 
A — 3.3, and the amplitudes are not directly compara- 
ble due to the different models and a slightly different 
definition of the amplitude. 

In conclusion, neither numerical errors nor shock for- 
mation can be the reason for the decay of the r mode, and 
its energy has to be transferred elsewhere, in particular 
into differential rotation. Wave-breaking does not occur 
in our models, but might play a role for even larger am- 
plitudes. 

Next, we try to get a more local picture of the r- 
mode decay. For this we study the integrand of the 
I = m = 2 current multipole. However, we first have to 
eliminate the influence of differential rotation, because 
locally the current multipole integrand has a significant 
(/)-component and therefore couples strongly to differen- 
tial rotation (and rotation as such when the background 
model is not subtracted), although this cancels out af- 
ter integration. It is therefore difficult to interpret plots 
based on the integrand itself, as done in [32]. Instead, we 
integrate along the ^-direction to get rid of any contribu- 
tion from perturbations with angular dependency m ^ 2. 
Snapshots at different times are shown in Fig. [12] While 
the total amplitude decays, the pattern is deformed, but 
not destroyed. This seems to contradict [32], where a 
breakdown of the mode pattern is reported. However, it 
is not obvious to what degree Fig. 5 in |32| is determined 
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FIG. 12: Time development of the (;6-integrated I — m = 
2 current multipole integrand for a nonlinearly decaying r 
mode of model MB85, with initial amplitude A — 3. (a)- 
(d) Absolute value at different times, (d) Complex phase 
difference to some reference point near the end of the decay. 
The phase was computed such that the usual phase-jump from 
— TT to TT is avoided, keeping the phase continuous. 



by the developing differential rotation instead of the r 
mode. 

We also found that the complex phase, which is spa- 
tially constant in the linear regime, develops a significant 
variance during the decay. Fig. [12] shows the variance at 
a late stage. This seems to imply that parts of the star 
oscillate out of phase. Note however that the local am- 
plitude goes to zero at the axis roughly quadratically, 
and the phase varies most strongly near the axis where it 
is very sensitive to contributions from secondary m — 2 
modes. We cannot rule out the possibility that the vari- 
ance is an artifact caused by other perturbations. It is 
worth mentioning that a large phase variance is reported 
in 'ST as well, although this was under the infiuence of an 
artificially large gravitational back-reaction force driving 
the r mode. 

In order to determine the importance of the phase vari- 
ance for our estimate of the r-mode energy from the mul- 
tipole moment, we recomputed the total multipole mo- 
ment, but used the absolute value of the complex inte- 
grand. The rationale is that the energy, in contrast to 
the multipole moment, is not sensitive to an axisymmet- 
ric phase shift. As shown in Fig. 
negligible. 



10 the differences are 



We note that the phase variance cannot be attributed 
to high amplitudes, since it is still present when the am- 
plitude has already decayed to values in the linear regime. 
If we assume for a moment that the phase variance of the 
current multipole really means that the oscillation of the 
velocity field is spatially out of phase, it follows that some 
other perturbation must be present, somehow infiuencing 
the r mode. It is worth noting that the differential ro- 



D. Search for mode coupling 



We now discuss possible nonlinear interactions of the r 
mode with secondary modes. For this, we looked at the 
Fourier spectra of the time evolution at various sample 
points. To study the time evolution of the spectra, we 
computed separate spectra for the first and second half 
of the evolution. Fig. 13 shows the spectra for the ve- 



locity components , for an initial amplitude A = 3. 
As one can see, the r mode is the dominant contribu- 
tion. However, there are significant peaks corresponding 
to oscillation modes in the inertial mode spectrum. We 
cannot identify those modes since the frequency resolu- 
tion of our spectra is insufficient to distinguish modes 
in the dense inertial mode spectrum. Interestingly, there 
are also peaks in the frequency range where the 2nd order 
scalar partial differential equation describing mode oscil- 
lations is of mixed elliptic-hyperbolic type, as discussed 
in [iO . It is unclear what the structure of solutions in 
this range would be, and if such solutions exist at all. 
It is however possible that the developing differential ro- 
tation shifts the location of this band. In any case, the 
only peaks we identified beside the r mode are the quasi- 
radial F mode and the axisymmetric /-mode, both quite 
insignificant for the velocity field. In the density, they are 
more visible, since the density perturbation of r modes 
in relation to the velocity perturbation is smaller than 
for pressure modes. 

None of the secondary modes with significant ampli- 
tude seem to grow. Only the increasing differential rota- 
tion is clearly visible in the spectra of v'^ (not plotted). 



From Fig. 13 we cannot confirm any mode coupling ef- 
fect, at least not of the magnitude found in [23l[32]. Note 
however that the spectra in [221 [31] are for an initial am- 
plitude A = 1.6, where the decay is slower. Looking 
at spectra from our simulations at lower amplitudes, we 
sometimes see growing peaks, but their amplitudes are 
tiny compared to the r mode. There are several differ- 
ent explanations for the discrepancies. First, the mode 
coupling for our models might saturate already during 
the timespan covered by the early stage Fourier spectra, 
thus being indistinguishable from secondary oscillations 
excited by the initial perturbation. Second, the mode 
coupling reported in [531 [S3 might not be the only cause 
of r-mode decay, and just happens to be less prominent 
for our models. Also, we might have chosen the wrong 
sample points, where an important secondary mode hap- 
pens to be small (Although we also studied some global 
quantities like multipole moments) . Our results thus can- 
not completely rule out mode coupling as the cause of the 
r-mode decay. 
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FIG. 13: Fourier spectra of the time evolution of « , u'" at 
a sample point x = y = z = 0.24 of model MB85, dur- 
ing the r-mode decay with initial amplitude ^ = 3. Shown 
are two spectra corresponding to the first and second half of 
the evolution. For comparison, vertical lines mark the known 
frequencies of the r mode, the quasi-radial F mode, and the 
axisymmetric / mode. The shaded area marks the transi- 
tional band between inertial and pressure modes (see main 
text). 



VI. SUMMARY AND DISCUSSION 

This work provides new evidence that r modes (with 
? = m = 2) of uniformly rotating neutron stars with a 
barotropic EOS decay rapidly if their dimensionless am- 
pHtude exceeds a model dependent value of order unity. 
The speed of decay depends only on the initial amplitude, 
and the decay does not stop even when the amplitude 
becomes small compared to unity. The r mode decays 
more slowly for higher rotation rates (at a fixed central 
density). Together with the decay of the r mode, strong 
differential rotation develops. The final rotation profile 
depends roughly on the distance to the axis. Close to the 
axis, the rotation is slowed down by a factor 0.5, while 
near the equator we observed speedups around 1.2. 

Our results are the first ones obtained using the rel- 
ativistic Cowling approximation, and are in good agree- 
ment with previous studies |311 132] of uniformly rotating 
stars which have been treated in the Newtonian frame- 
work, but without artificially fixing the gravitational 
field. Thus, the cause is unrelated both to relativistic ef- 
fects and to the changes of the gravitational field induced 
by the fluid. Those studies also used a different way to ex- 
cite the r mode. While they used an artificially increased 
gravitational radiation reaction force to slowly drive the 
r mode to large amplitudes, we perturbed the initial data 
with linearly scaled exact eigenfunctions. This is note- 
worthy because it implies that the decay and the dif- 
ferential rotation are not sensitive to the amount and 
composition of other modes in the initial data. 

However, we also found some differences. The afore- 
mentioned Newtonian studies reported the occurrence of 
either wave breaking or strong mode coupling together 



with the decay. We found no wave-breaking, although the 
oscillations near the surface are definitely in the nonlinear 
regime. This is not a contradiction. Given the different 
models, our maximum amplitude, although comparable 
to the the reported wave-breaking case, was probably just 
not large enough. Nevertheless, wave-breaking is not nec- 
essary for the r-mode decay. We also cannot confirm the 
presence of significant mode coupling. However, since 
we only analyzed the time evolution of selected sample 
points and a few multipole moments, we cannot rule it 
out completely. 

We also studied the energy budget of the process. For 
this, we derived a measure for the energy of fluid modes, 
which estimates the energy difference between a state 
perturbed with the linear eigenfunction and a ground 
state with the same angular momentum profile and total 
mass. For this we took advantage of the fact that us- 
ing an artificially fixed axisymmetric spacetime implies 
the existence of a conserved energy and angular momen- 
tum besides the conserved mass. During the nonlinear 
decay, we observed a significant loss of conserved energy 
caused mostly, if not completely, by numerical errors, as 
we found from convergence tests. The energy loss of the 
r mode however was greater, and not sensitive to the 
numerical resolution. We can thus conclude that the en- 
ergy of the r mode is not dissipated directly, e.g. due 
to wave breaking, shock formation, or numerical errors. 
Instead it is converted mostly into differential rotation, 
which also causes a deformation of the star, increasing 
the equatorial radius by a few percent. 

Last but not least, we provided eigenfunctions and fre- 
quencies of r modes in the relativistic Cowling approx- 
imation for rapidly rotating stars. We found excellent 
agreement with the frequencies found in [16 , thus val- 
idating those results. Within numerical accuracy, the 
eigenfunctions are smooth. Eigenfunctions and frequen- 
cies are still very similar to the Newtonian slow rota- 
tion case. We also estimated the gravitational luminos- 
ity, wave strain, and angular momentum loss caused by 
r modes, and found that the current multipole is still 
strongly dominant for the rapidly rotating case. For our 
models, an r mode with unit dimensionless amplitude at 
a distance of 10 Mpc causes a wave strain of the order 
10-24. 

Finally we would like to speculate a bit. Although the 
decay timescale diverges quickly with decreasing ampli- 
tude, we have no proof that the effect vanishes completely 
below some critical amplitude. Let us assume that the 
r-mode growth due to the CFS instability and the de- 
cay described in this work are balanced at some ampli- 
tude too small to be relevant as a source for detectable 
gravitational radiation. It is plausible to assume that 
the r-mode energy loss is still converted to differential 
rotation. The effect would then be cumulative as long 
as the CFS instability is active. It might well be that 
the CFS instability of the r mode does not induce large 
amplitude oscillations, but differential rotation of similar 
energy. This effect, which is of course purely hypothet- 
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ical, might be important for the time evolution of the 
rotation profile of newborn neutron stars, and also may 
lead to strong amplification of the magnetic field. 
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